\section*{Boundary Integral Method for Laplacian}

We wanna solve a 2D Laplacian problem with homogeneous Dirichlet boundary condition, which can be stated as follow

\begin{equation}\nonumber
\begin{cases}
\nabla^2 u(\textbf{x}) = 0, \textbf{x} \in \Omega\\
u(\textbf{x}) = g(\textbf{x}), \textbf{x} \in \partial \Omega
\end{cases}
\end{equation}

The way we wanna solve the above problem is by using Indirect Formulation of Boundary Integral Method, and we assume the final solution should look like a double layer potential with moment or charge density h, defined as follow

\begin{equation}\nonumber
\hat{u}(\textbf{x}) = -\oint_{\partial\Omega}h(\textbf{y})\frac{\partial\Phi}{\partial\textbf{n}_y}(\textbf{x} - \textbf{y})ds(\textbf{y})
\end{equation}

Here $\Phi$ is the fundamental solution to the Laplacian equation and h is some smooth function defined on the boundary $\partial \Omega$. Such function satisfies the following property automatically

\begin{equation}\nonumber
\nabla^2 \hat{u}(\textbf{x}) = 0, \textbf{x} \in \bar{\Omega}
\end{equation}

given a Laplacian problem with homogeneous Dirichlet boundary condition, we are looking for a solution by making use of the doule layer potential representation. Thus we can solve for the unknown charge density function h by enforcing it to satisfy the boundary constraints.

The field generated by the double layer potential is actually harmonic in the interior domain $\Omega$ and exterior domain $\Omega^{C}$, but has a jumping condition on the layer boundary. The limit from the interior domain and exterior domain are related to the jumping region as follow

\begin{equation}\nonumber
\begin{cases}
\displaystyle{\lim_{\textbf{x} \to \textbf{x}_0}} \hat{u}(\textbf{x}) = \frac{1}{2}h(\textbf{x}_0) + \hat{u}(\textbf{x}_0), \textbf{x} \in \bar{\Omega}\\
\displaystyle{\lim_{\textbf{x} \to \textbf{x}_0}} \hat{u}(\textbf{x}) = -\frac{1}{2}h(\textbf{x}_0) + \hat{u}(\textbf{x}_0), \textbf{x} \in \Omega^{C}
\end{cases}
\end{equation}

Then to solve our original problem, we set

\begin{equation}\nonumber
g(\textbf{x}_0) = \displaystyle{\lim_{\textbf{x} \to \textbf{x}_0}} \hat{u}(\textbf{x}) = \frac{1}{2}h(\textbf{x}_0) + \hat{u}(\textbf{x}_0), \textbf{x} \in \bar{\Omega}
\end{equation}

to resolve the charge density function $h(\cdot)$. Here we also list the analytical representation for the fundamental solution and its directional derivative

\begin{equation}\nonumber
\begin{split}
\Phi(\textbf{x} - \textbf{y}) = 
\begin{cases}
-\frac{1}{2\pi}\ln|\textbf{x} - \textbf{y}| \quad\quad\quad\quad\quad n=2\\
\frac{1}{n(n-2)\alpha(n)}\cdot\frac{1}{|\textbf{x}-\textbf{y}|^{n-2}} \,\,\,\quad\quad n\geq 3
\end{cases}\\
\frac{\partial\Phi}{\partial \nu_{\textbf{y}}}(\textbf{x} - \textbf{y}) = \nabla_{\textbf{y}}\Phi(\textbf{x} - \textbf{y})\cdot \nu(\textbf{y}) = \frac{(\textbf{x} - \textbf{y})\cdot\nu(\textbf{y})}{n\alpha(n)|\textbf{x} - \textbf{y}|}
\end{split}
\end{equation}

where $\alpha(n)$ is the volume of unit sphere in $\mathbb{R}^n$ and $\nu$ is a direction vector represented in the frame of \textbf{y}.

In the discrete setting, one can apply Galerkin scheme. Here we choose a collocated scheme, in which case we divide the boundary $\partial \Omega$ into nonoverlapping elements like line segments in 2D or triangles in 3D. From an approximation point of view, we define an approximate solution space  on the boundary $\aleph^{\infty}(\partial \Omega)$ made of several basis functions, defined as follow

\begin{equation}\nonumber
\aleph^{\infty}(\partial\Omega) = \{\Psi(\textbf{x}) | \Psi(\textbf{x}) \in L_{\infty}(\partial \Omega), \textbf{x} \in \partial \Omega\}
\end{equation}

Suppose we pick a set of basis functions just from $\aleph^{\infty}(\partial\Omega)$ and they span a subspace of $\aleph^{\infty}(\partial\Omega)$, let's call it $\aleph^{\infty}_{H}(\partial\Omega)$ . The the true charge density function can be approximated by a weighted sum of these basis functions

\begin{equation}\nonumber
h(\textbf{x}) = \tilde{h}(\textbf{x}) = \sum_j h_j \Psi_j(\textbf{x})
\end{equation}

Suppose these basis functions $\Psi_j$ are nodal basis functions and we evaluate the boundary value at those nodal points, which will give us

\begin{equation}\nonumber
\begin{split}
g(\textbf{x}_i) &= \frac{1}{2}h(\textbf{x}_i) - \oint_{\partial\Omega}\sum_{j}h_j\Psi_j(\textbf{y})\frac{\partial\Phi}{\partial\textbf{n}_y}(\textbf{x}_i - \textbf{y})ds(\textbf{y})\\
&= \frac{1}{2}h(\textbf{x}_i) -  \sum_{j}h_j [\oint_{\partial\Omega}\Psi_j(\textbf{y})\frac{\partial\Phi}{\partial\textbf{n}_y}(\textbf{x}_i - \textbf{y})ds(\textbf{y})]
\end{split}
\end{equation}

This leads to a linear system and if we pick enough points, we will resolve the weights for each basis function.

\begin{equation}\nonumber
\textbf{Ah} = \textbf{g}, \quad\quad \textbf{A}_{ij} = 
\begin{cases}
-\displaystyle{\oint_{\partial\Omega}\Psi_j(\textbf{y})\frac{\partial\Phi}{\partial\textbf{n}_y}(\textbf{x}_i - \textbf{y})ds(\textbf{y})}, \quad\quad\quad\quad\quad i \neq j\\
\frac{1}{2} - \displaystyle{\oint_{\partial\Omega}\Psi_j(\textbf{y})\frac{\partial\Phi}{\partial\textbf{n}_y}(\textbf{x}_i - \textbf{y})ds(\textbf{y})}, \quad\quad\quad\quad i = j
\end{cases}
\end{equation}